home *** CD-ROM | disk | FTP | other *** search
/ SGI Developer Toolbox 6.1 / SGI Developer Toolbox 6.1 - Disc 4.iso / lib / mathlib / libblas / src_original / dsbmv.f < prev    next >
Encoding:
Text File  |  1994-08-02  |  9.7 KB  |  307 lines
  1. *
  2. ************************************************************************
  3. *
  4.       SUBROUTINE DSBMV ( UPLO, N, K, ALPHA, A, LDA, X, INCX,
  5.      $                   BETA, Y, INCY )
  6. *     .. Scalar Arguments ..
  7.       DOUBLE PRECISION   ALPHA, BETA
  8.       INTEGER            INCX, INCY, K, LDA, N
  9.       CHARACTER*1        UPLO
  10. *     .. Array Arguments ..
  11.       DOUBLE PRECISION   A( LDA, * ), X( * ), Y( * )
  12. *     ..
  13. *
  14. *  Purpose
  15. *  =======
  16. *
  17. *  DSBMV  performs the matrix-vector  operation
  18. *
  19. *     y := alpha*A*x + beta*y,
  20. *
  21. *  where alpha and beta are scalars, x and y are n element vectors and
  22. *  A is an n by n symmetric band matrix, with k super-diagonals.
  23. *
  24. *  Parameters
  25. *  ==========
  26. *
  27. *  UPLO   - CHARACTER*1.
  28. *           On entry, UPLO specifies whether the upper or lower
  29. *           triangular part of the band matrix A is being supplied as
  30. *           follows:
  31. *
  32. *              UPLO = 'U' or 'u'   The upper triangular part of A is
  33. *                                  being supplied.
  34. *
  35. *              UPLO = 'L' or 'l'   The lower triangular part of A is
  36. *                                  being supplied.
  37. *
  38. *           Unchanged on exit.
  39. *
  40. *  N      - INTEGER.
  41. *           On entry, N specifies the order of the matrix A.
  42. *           N must be at least zero.
  43. *           Unchanged on exit.
  44. *
  45. *  K      - INTEGER.
  46. *           On entry, K specifies the number of super-diagonals of the
  47. *           matrix A. K must satisfy  0 .le. K.
  48. *           Unchanged on exit.
  49. *
  50. *  ALPHA  - DOUBLE PRECISION.
  51. *           On entry, ALPHA specifies the scalar alpha.
  52. *           Unchanged on exit.
  53. *
  54. *  A      - DOUBLE PRECISION array of DIMENSION ( LDA, n ).
  55. *           Before entry with UPLO = 'U' or 'u', the leading ( k + 1 )
  56. *           by n part of the array A must contain the upper triangular
  57. *           band part of the symmetric matrix, supplied column by
  58. *           column, with the leading diagonal of the matrix in row
  59. *           ( k + 1 ) of the array, the first super-diagonal starting at
  60. *           position 2 in row k, and so on. The top left k by k triangle
  61. *           of the array A is not referenced.
  62. *           The following program segment will transfer the upper
  63. *           triangular part of a symmetric band matrix from conventional
  64. *           full matrix storage to band storage:
  65. *
  66. *                 DO 20, J = 1, N
  67. *                    M = K + 1 - J
  68. *                    DO 10, I = MAX( 1, J - K ), J
  69. *                       A( M + I, J ) = matrix( I, J )
  70. *              10    CONTINUE
  71. *              20 CONTINUE
  72. *
  73. *           Before entry with UPLO = 'L' or 'l', the leading ( k + 1 )
  74. *           by n part of the array A must contain the lower triangular
  75. *           band part of the symmetric matrix, supplied column by
  76. *           column, with the leading diagonal of the matrix in row 1 of
  77. *           the array, the first sub-diagonal starting at position 1 in
  78. *           row 2, and so on. The bottom right k by k triangle of the
  79. *           array A is not referenced.
  80. *           The following program segment will transfer the lower
  81. *           triangular part of a symmetric band matrix from conventional
  82. *           full matrix storage to band storage:
  83. *
  84. *                 DO 20, J = 1, N
  85. *                    M = 1 - J
  86. *                    DO 10, I = J, MIN( N, J + K )
  87. *                       A( M + I, J ) = matrix( I, J )
  88. *              10    CONTINUE
  89. *              20 CONTINUE
  90. *
  91. *           Unchanged on exit.
  92. *
  93. *  LDA    - INTEGER.
  94. *           On entry, LDA specifies the first dimension of A as declared
  95. *           in the calling (sub) program. LDA must be at least
  96. *           ( k + 1 ).
  97. *           Unchanged on exit.
  98. *
  99. *  X      - DOUBLE PRECISION array of DIMENSION at least
  100. *           ( 1 + ( n - 1 )*abs( INCX ) ).
  101. *           Before entry, the incremented array X must contain the
  102. *           vector x.
  103. *           Unchanged on exit.
  104. *
  105. *  INCX   - INTEGER.
  106. *           On entry, INCX specifies the increment for the elements of
  107. *           X. INCX must not be zero.
  108. *           Unchanged on exit.
  109. *
  110. *  BETA   - DOUBLE PRECISION.
  111. *           On entry, BETA specifies the scalar beta.
  112. *           Unchanged on exit.
  113. *
  114. *  Y      - DOUBLE PRECISION array of DIMENSION at least
  115. *           ( 1 + ( n - 1 )*abs( INCY ) ).
  116. *           Before entry, the incremented array Y must contain the
  117. *           vector y. On exit, Y is overwritten by the updated vector y.
  118. *
  119. *  INCY   - INTEGER.
  120. *           On entry, INCY specifies the increment for the elements of
  121. *           Y. INCY must not be zero.
  122. *           Unchanged on exit.
  123. *
  124. *
  125. *  Level 2 Blas routine.
  126. *
  127. *  -- Written on 22-October-1986.
  128. *     Jack Dongarra, Argonne National Lab.
  129. *     Jeremy Du Croz, Nag Central Office.
  130. *     Sven Hammarling, Nag Central Office.
  131. *     Richard Hanson, Sandia National Labs.
  132. *
  133. *
  134. *     .. Parameters ..
  135.       DOUBLE PRECISION   ONE         , ZERO
  136.       PARAMETER        ( ONE = 1.0D+0, ZERO = 0.0D+0 )
  137. *     .. Local Scalars ..
  138.       DOUBLE PRECISION   TEMP1, TEMP2
  139.       INTEGER            I, INFO, IX, IY, J, JX, JY, KPLUS1, KX, KY, L
  140. *     .. External Functions ..
  141.       LOGICAL            LSAME
  142.       EXTERNAL           LSAME
  143. *     .. External Subroutines ..
  144.       EXTERNAL           XERBLA
  145. *     .. Intrinsic Functions ..
  146.       INTRINSIC          MAX, MIN
  147. *     ..
  148. *     .. Executable Statements ..
  149. *
  150. *     Test the input parameters.
  151. *
  152.       INFO = 0
  153.       IF     ( .NOT.LSAME( UPLO, 'U' ).AND.
  154.      $         .NOT.LSAME( UPLO, 'L' )      )THEN
  155.          INFO = 1
  156.       ELSE IF( N.LT.0 )THEN
  157.          INFO = 2
  158.       ELSE IF( K.LT.0 )THEN
  159.          INFO = 3
  160.       ELSE IF( LDA.LT.( K + 1 ) )THEN
  161.          INFO = 6
  162.       ELSE IF( INCX.EQ.0 )THEN
  163.          INFO = 8
  164.       ELSE IF( INCY.EQ.0 )THEN
  165.          INFO = 11
  166.       END IF
  167.       IF( INFO.NE.0 )THEN
  168.          CALL XERBLA( 'DSBMV ', INFO )
  169.          RETURN
  170.       END IF
  171. *
  172. *     Quick return if possible.
  173. *
  174.       IF( ( N.EQ.0 ).OR.( ( ALPHA.EQ.ZERO ).AND.( BETA.EQ.ONE ) ) )
  175.      $   RETURN
  176. *
  177. *     Set up the start points in  X  and  Y.
  178. *
  179.       IF( INCX.GT.0 )THEN
  180.          KX = 1
  181.       ELSE
  182.          KX = 1 - ( N - 1 )*INCX
  183.       END IF
  184.       IF( INCY.GT.0 )THEN
  185.          KY = 1
  186.       ELSE
  187.          KY = 1 - ( N - 1 )*INCY
  188.       END IF
  189. *
  190. *     Start the operations. In this version the elements of the array A
  191. *     are accessed sequentially with one pass through A.
  192. *
  193. *     First form  y := beta*y.
  194. *
  195.       IF( BETA.NE.ONE )THEN
  196.          IF( INCY.EQ.1 )THEN
  197.             IF( BETA.EQ.ZERO )THEN
  198.                DO 10, I = 1, N
  199.                   Y( I ) = ZERO
  200.    10          CONTINUE
  201.             ELSE
  202.                DO 20, I = 1, N
  203.                   Y( I ) = BETA*Y( I )
  204.    20          CONTINUE
  205.             END IF
  206.          ELSE
  207.             IY = KY
  208.             IF( BETA.EQ.ZERO )THEN
  209.                DO 30, I = 1, N
  210.                   Y( IY ) = ZERO
  211.                   IY      = IY   + INCY
  212.    30          CONTINUE
  213.             ELSE
  214.                DO 40, I = 1, N
  215.                   Y( IY ) = BETA*Y( IY )
  216.                   IY      = IY           + INCY
  217.    40          CONTINUE
  218.             END IF
  219.          END IF
  220.       END IF
  221.       IF( ALPHA.EQ.ZERO )
  222.      $   RETURN
  223.       IF( LSAME( UPLO, 'U' ) )THEN
  224. *
  225. *        Form  y  when upper triangle of A is stored.
  226. *
  227.          KPLUS1 = K + 1
  228.          IF( ( INCX.EQ.1 ).AND.( INCY.EQ.1 ) )THEN
  229.             DO 60, J = 1, N
  230.                TEMP1 = ALPHA*X( J )
  231.                TEMP2 = ZERO
  232.                L     = KPLUS1 - J
  233.                DO 50, I = MAX( 1, J - K ), J - 1
  234.                   Y( I ) = Y( I ) + TEMP1*A( L + I, J )
  235.                   TEMP2  = TEMP2  + A( L + I, J )*X( I )
  236.    50          CONTINUE
  237.                Y( J ) = Y( J ) + TEMP1*A( KPLUS1, J ) + ALPHA*TEMP2
  238.    60       CONTINUE
  239.          ELSE
  240.             JX = KX
  241.             JY = KY
  242.             DO 80, J = 1, N
  243.                TEMP1 = ALPHA*X( JX )
  244.                TEMP2 = ZERO
  245.                IX    = KX
  246.                IY    = KY
  247.                L     = KPLUS1 - J
  248.                DO 70, I = MAX( 1, J - K ), J - 1
  249.                   Y( IY ) = Y( IY ) + TEMP1*A( L + I, J )
  250.                   TEMP2   = TEMP2   + A( L + I, J )*X( IX )
  251.                   IX      = IX      + INCX
  252.                   IY      = IY      + INCY
  253.    70          CONTINUE
  254.                Y( JY ) = Y( JY ) + TEMP1*A( KPLUS1, J ) + ALPHA*TEMP2
  255.                JX      = JX      + INCX
  256.                JY      = JY      + INCY
  257.                IF( J.GT.K )THEN
  258.                   KX = KX + INCX
  259.                   KY = KY + INCY
  260.                END IF
  261.    80       CONTINUE
  262.          END IF
  263.       ELSE
  264. *
  265. *        Form  y  when lower triangle of A is stored.
  266. *
  267.          IF( ( INCX.EQ.1 ).AND.( INCY.EQ.1 ) )THEN
  268.             DO 100, J = 1, N
  269.                TEMP1  = ALPHA*X( J )
  270.                TEMP2  = ZERO
  271.                Y( J ) = Y( J )       + TEMP1*A( 1, J )
  272.                L      = 1            - J
  273.                DO 90, I = J + 1, MIN( N, J + K )
  274.                   Y( I ) = Y( I ) + TEMP1*A( L + I, J )
  275.                   TEMP2  = TEMP2  + A( L + I, J )*X( I )
  276.    90          CONTINUE
  277.                Y( J ) = Y( J ) + ALPHA*TEMP2
  278.   100       CONTINUE
  279.          ELSE
  280.             JX = KX
  281.             JY = KY
  282.             DO 120, J = 1, N
  283.                TEMP1   = ALPHA*X( JX )
  284.                TEMP2   = ZERO
  285.                Y( JY ) = Y( JY )       + TEMP1*A( 1, J )
  286.                L       = 1             - J
  287.                IX      = JX
  288.                IY      = JY
  289.                DO 110, I = J + 1, MIN( N, J + K )
  290.                   IX      = IX      + INCX
  291.                   IY      = IY      + INCY
  292.                   Y( IY ) = Y( IY ) + TEMP1*A( L + I, J )
  293.                   TEMP2   = TEMP2   + A( L + I, J )*X( IX )
  294.   110          CONTINUE
  295.                Y( JY ) = Y( JY ) + ALPHA*TEMP2
  296.                JX      = JX      + INCX
  297.                JY      = JY      + INCY
  298.   120       CONTINUE
  299.          END IF
  300.       END IF
  301. *
  302.       RETURN
  303. *
  304. *     End of DSBMV .
  305. *
  306.       END
  307.